Metagenomics Workshop, November 9, 2022

Introduction and Background

Agenda

Part I - Setting up the Data

• Define the differences between Metabarcoding/Amplicon sequencing and Whole Shotgun Sequencing Metagenomics
• Get some practice in RStudio OnDemand and using the Tufts Cluster
• Learn about accessing interesting NCBI data through BioProjects and SRA Run Selector
• Refine our quality control best practices to preprocess data we have generated or downloaded
• Fastqc
• Multiqc
• Kmer-based filtering and trimming with BCBio tools

Part II - Practice a typical metagenomic workflow

• Metagenome Assembly with Megahit
• Assessing the Quality of an Assembly with Quast
• Preliminary contamination screening with kraken
• Visualization of kraken reports with the krona tool

The following items are not hands on, but will be shown as part of use case of a recent Metagenomics Workflow done in collaboration with a Tufts Researcher

• Metagenome Binning and Refinement with MetaBat2, MaxBin and Concoct
• Visualization of Sequence of Metagenomic Bins with BlobTools
• Extracting Relevant Reads with BWA
• Reassembly with Megahit from the Bins
• Taxonomic Classification with Kaiju
• Functional Classification with Prokka

If this sounds like a lot of ground to cover, you are correct!

In some portions of this tutorial, we don’t actually run the code, but jump ahead to the next step.

For these “skipped” portions, the code to run the process is provided.

---

Attach public library paths to your RStudio OnDemand session

The setup code adds the libraries in the order they will be searced for packages.

Setting .libpaths() to Tufts shared folders can reduce the amount of time spent downloading packages.

Please reach out to us about any difficult or tricky installs, sometimes it is a library package error that has to be fixed by an TTS HPC Research Technology Specialist.

Please make sure that the output to this code chunk has the shared library, “/cluster/tufts/hpc/tools/R/4.0.0” in position [1]

Example output (it depends what else you have been doing in RStudio before today).

[1] “/cluster/tufts/hpc/tools/R/4.0.0” [2] “/cluster/home/arhode05/R/x86_64-pc-linux-gnu-library/4.0” [3] “/opt/shared/R/4.0.0/lib64/R/library”

Commands to specify the shared directory as the first place to look for R libraries and check the output.

.libPaths(c('/cluster/tufts/hpc/tools/R/4.0.0',.libPaths()))

.libPaths()

This Metagenomics Tutorial is mainly a set of Bash Commands

This tutorial uses code chunks. To run an individual code chunk while the notebook is open in OnDemand RStudio, please press the green triangle on the upper right corner of the code chunk.

You may be wondering why we went through the trouble of having bash commands inside an R notebook for this tutorial. The simplest answer is that RStudio has the handy ability to open visualization files with a simple click on the file name in the lower right screen. This saves a lot of time and confusing navigation when doing this type of project, that produces a lot of intermediate files.

RStudio can run code chunks in other languages by resetting the code chunks to read that language.

• {bash} indicates that the code chunk is in command line syntax, so these commands can be copied and pasted directly into a terminal

• {r} indicates that the code chunk uses R-code

To save a bunch of copying/pasting, we are just going to run both types of these chunks inside our notebook using RStudio OnDemand.

The order of chunks is important, so if you skip one or get lost, it is possible to use the dropdown menu next to the word “Run” at the top of this window to run multiple chunks to get you back to where you started.

Any questions?

---

What other “kernels” are available in R notebook code chunks? (Hint - click the down arrow next to the "Insert Chunk command at the top of this window)

---

Modules and Conda Environments

In this workshop, we will sometimes be loading modules from the HPC to use in a code chunk.

Here is an example code chunk that loads a module for samtools and runs the command with the help flag to generate some output to test that the module loaded correctly.

If the bash environment is not designated in the setup chunk, you may need to add this bash command at the beginning of the chunk to enable module loading.

source ~/.bashrc

module load samtools/1.9
samtools --help

The behavior of modules in the code chunks differs from our direct command line.

The Rmarkdown notebook will not remember it from chunk to chunk. So load the modules in the chunks where the commands are needed.

An example with conda environments will be shown in Part II.

---

What happens when you run this chunk?

module list

No Modulefiles Currently Loaded.

---

Why use modules instead of conda environments?

• Conda environments are great for installing specific versions of programs, but often bring in clashing dependencies, and only one is active at a time.

• Using HPC modules can provide more integration because you can load several at one time.

• Of course, conda environments could be used if you load all the required programs into the same environment. We will see an example of this later with Metawrap

Introducing the Data - Giant Panda ticks

Ticks are capable of spreading pathogens (viruses, bacteria and other parasites) among a host population. Not much is known about the spread of tick-borne viruses in Giant Pandas. The density of giant pandas in breeding captivity in Sichuan Province, China may lead to more transmission of these diseases. The researchers who loaded this data into the SRA were interested in characterizing the virome of these ticks and to compare these with the virome of the pandas.

We are not going to restrict our analysis to the virome, and just take a general look at what these ticks are carryng around in general using practical metagenomics workflows.

Check out the bioproject here

Check out the final paper here

Is he sleepy or is he sick?

We are going to bring in some of the data from this bioproject to our Metagenomics Tutorial space.

What information can we find in the metadata?

SRA Metadata - every project is different!

Example from BioProject Microbes From Mum

Mum sponge metadata

The Metadata from our BioProject

Panda tick metadata

Question - Why do most of the reads of “spot length” equal to 502?

Something you may have noticed is that the reads are suspicially even in number, at length 502. The experimental design states that the reads were paired-end 2x250.

Question - Who is responsible for loading BioProject metadata to the NCBI?

Obtaining Data from SRA

The data we are going to copy over can be obtained directly from the SRA by using our command line module for the SRAToolkit. This can take a little while, depending on the internet connection and which source the files are pulled from (on-prem, AWS or GCP).

“Spots” in SRA are not the same as “reads”

• A spot is all the info you got from one “spot” on the flow cell.
• You get 4 reads per spot with today’s illumina sequencing: forward barcode, forward read, reverse barcode, reverse read.

NCBI SRA “Spots”

Practice downloading data from SRA

We have a toolkit on the HPC cluster that allows users to pull data directly from SRA. The first step is configuring the SRA toolkit with a unique user ID. The second step is running either fastq-dump or fasterq-dump. We will use fastq-dump to demonstrate the options available.

module load /cluster/tufts/bio/tools/module_files/sra/3.0.0

#The next few lines set up your configuration for SRA and is necessary if you have not used the SRA toolkit previously.

#The vdb-config is an interactive tool, to generate a unique identifier. Since it is difficult to use it interactively from ROnDemand, this command opens the interactive screen for three seconds then closes it, which generates the user id.

vdb-config -i & read -t 3 ; kill $!

[?25l[?1000h[?1002h

module load /cluster/tufts/bio/tools/module_files/sra/3.0.0

mkdir -p test_SRA
fastq-dump --outdir test_SRA -X 500 --skip-technical --read-filter pass --readids --minReadLen 200 --split-e SRR18086364

Rejected 149 READS because READLEN < 200
Read 500 spots for SRR18086364
Written 426 spots for SRR18086364

Command Line Parameters

• Grap first million lines -X 500
• Only “biological reads” --skip-technical
• Remove reads that are mostly N’s --read-filter pass
• Append .1 to forward read header .2 to reverse read header –readids
• Keep mated pairs and singletons, use -–split-e and a three files will be generated if orphans are present (forward, reverse, and unmatched mates)

There should be three new files in your Metagenomics2022/test_SRA directory.

The output is small enough that we can look into the fastq file by clicking on it.

• Confirm the .1 and .2 were appended to your forward read headers and reverse read headers
• What happened inside the singleton file? Do you see .1, .2 or both?
• QC Best practice: Count the number of lines in the output files

wc -l ./test_SRA/*.fastq

  1700 ./test_SRA/SRR18086364_pass_1.fastq
  1700 ./test_SRA/SRR18086364_pass_2.fastq
     4 ./test_SRA/SRR18086364_pass.fastq
  3404 total

• How many lines are in our fastq output files and why?

  (Hint 1: we asked for 500 “spots”, not reads) (Hint 2: how did our parameter choices affect our outputs?)

A Note about Fasterq-dump

In some cases, the newer version of this tool fasterq-dump will automatically run the split-e version of the command without adding this parameter. It is worth looking at the options before using it, because some of the parameter meanings have changed.

The following codeblock is for reference, to use the chunk, remove the # at the beginning of the line.

#/cluster/tufts/hpc/tools/spack/linux-rhel7-ivybridge/gcc-9.3.0/fastqc-0.11.9-eeij5jwnwau3ttx4m53bxafmxbhm4fjz/bin/fastqc ./MGData/*overrep.fastq

#source ~/.bashrc
module load fastqc
mkdir -p fastqc
fastqc ./data/*pass*.fastq -o fastqc

Started analysis of SRR18085691_pass_1.fastq
Approx 5% complete for SRR18085691_pass_1.fastq
Approx 10% complete for SRR18085691_pass_1.fastq
Approx 15% complete for SRR18085691_pass_1.fastq
Approx 20% complete for SRR18085691_pass_1.fastq
Approx 25% complete for SRR18085691_pass_1.fastq
Approx 30% complete for SRR18085691_pass_1.fastq
Approx 35% complete for SRR18085691_pass_1.fastq
Approx 40% complete for SRR18085691_pass_1.fastq
Approx 45% complete for SRR18085691_pass_1.fastq
Approx 50% complete for SRR18085691_pass_1.fastq
Approx 55% complete for SRR18085691_pass_1.fastq
Approx 60% complete for SRR18085691_pass_1.fastq
Approx 65% complete for SRR18085691_pass_1.fastq
Approx 70% complete for SRR18085691_pass_1.fastq
Approx 75% complete for SRR18085691_pass_1.fastq
Approx 80% complete for SRR18085691_pass_1.fastq
Approx 85% complete for SRR18085691_pass_1.fastq
Approx 90% complete for SRR18085691_pass_1.fastq
Approx 95% complete for SRR18085691_pass_1.fastq

Analysis complete for SRR18085691_pass_1.fastq

Approx 100% complete for SRR18085691_pass_1.fastq
Started analysis of SRR18085691_pass_2.fastq
Approx 5% complete for SRR18085691_pass_2.fastq
Approx 10% complete for SRR18085691_pass_2.fastq
Approx 15% complete for SRR18085691_pass_2.fastq
Approx 20% complete for SRR18085691_pass_2.fastq
Approx 25% complete for SRR18085691_pass_2.fastq
Approx 30% complete for SRR18085691_pass_2.fastq
Approx 35% complete for SRR18085691_pass_2.fastq
Approx 40% complete for SRR18085691_pass_2.fastq
Approx 45% complete for SRR18085691_pass_2.fastq
Approx 50% complete for SRR18085691_pass_2.fastq
Approx 55% complete for SRR18085691_pass_2.fastq
Approx 60% complete for SRR18085691_pass_2.fastq
Approx 65% complete for SRR18085691_pass_2.fastq
Approx 70% complete for SRR18085691_pass_2.fastq
Approx 75% complete for SRR18085691_pass_2.fastq
Approx 80% complete for SRR18085691_pass_2.fastq
Approx 85% complete for SRR18085691_pass_2.fastq
Approx 90% complete for SRR18085691_pass_2.fastq
Approx 95% complete for SRR18085691_pass_2.fastq

Analysis complete for SRR18085691_pass_2.fastq

Approx 100% complete for SRR18085691_pass_2.fastq
Started analysis of SRR18085708_pass_1.fastq
Approx 5% complete for SRR18085708_pass_1.fastq
Approx 10% complete for SRR18085708_pass_1.fastq
Approx 15% complete for SRR18085708_pass_1.fastq
Approx 20% complete for SRR18085708_pass_1.fastq
Approx 25% complete for SRR18085708_pass_1.fastq
Approx 30% complete for SRR18085708_pass_1.fastq
Approx 35% complete for SRR18085708_pass_1.fastq
Approx 40% complete for SRR18085708_pass_1.fastq
Approx 45% complete for SRR18085708_pass_1.fastq
Approx 50% complete for SRR18085708_pass_1.fastq
Approx 55% complete for SRR18085708_pass_1.fastq
Approx 60% complete for SRR18085708_pass_1.fastq
Approx 65% complete for SRR18085708_pass_1.fastq
Approx 70% complete for SRR18085708_pass_1.fastq
Approx 75% complete for SRR18085708_pass_1.fastq
Approx 80% complete for SRR18085708_pass_1.fastq
Approx 85% complete for SRR18085708_pass_1.fastq
Approx 90% complete for SRR18085708_pass_1.fastq
Approx 95% complete for SRR18085708_pass_1.fastq

Analysis complete for SRR18085708_pass_1.fastq

Approx 100% complete for SRR18085708_pass_1.fastq
Started analysis of SRR18085708_pass_2.fastq
Approx 5% complete for SRR18085708_pass_2.fastq
Approx 10% complete for SRR18085708_pass_2.fastq
Approx 15% complete for SRR18085708_pass_2.fastq
Approx 20% complete for SRR18085708_pass_2.fastq
Approx 25% complete for SRR18085708_pass_2.fastq
Approx 30% complete for SRR18085708_pass_2.fastq
Approx 35% complete for SRR18085708_pass_2.fastq
Approx 40% complete for SRR18085708_pass_2.fastq
Approx 45% complete for SRR18085708_pass_2.fastq
Approx 50% complete for SRR18085708_pass_2.fastq
Approx 55% complete for SRR18085708_pass_2.fastq
Approx 60% complete for SRR18085708_pass_2.fastq
Approx 65% complete for SRR18085708_pass_2.fastq
Approx 70% complete for SRR18085708_pass_2.fastq
Approx 75% complete for SRR18085708_pass_2.fastq
Approx 80% complete for SRR18085708_pass_2.fastq
Approx 85% complete for SRR18085708_pass_2.fastq
Approx 90% complete for SRR18085708_pass_2.fastq
Approx 95% complete for SRR18085708_pass_2.fastq

Analysis complete for SRR18085708_pass_2.fastq

Approx 100% complete for SRR18085708_pass_2.fastq
Started analysis of SRR18086277_pass_1.fastq
Approx 5% complete for SRR18086277_pass_1.fastq
Approx 10% complete for SRR18086277_pass_1.fastq
Approx 15% complete for SRR18086277_pass_1.fastq
Approx 20% complete for SRR18086277_pass_1.fastq
Approx 25% complete for SRR18086277_pass_1.fastq
Approx 30% complete for SRR18086277_pass_1.fastq
Approx 35% complete for SRR18086277_pass_1.fastq
Approx 40% complete for SRR18086277_pass_1.fastq
Approx 45% complete for SRR18086277_pass_1.fastq
Approx 50% complete for SRR18086277_pass_1.fastq
Approx 55% complete for SRR18086277_pass_1.fastq
Approx 60% complete for SRR18086277_pass_1.fastq
Approx 65% complete for SRR18086277_pass_1.fastq
Approx 70% complete for SRR18086277_pass_1.fastq
Approx 75% complete for SRR18086277_pass_1.fastq
Approx 80% complete for SRR18086277_pass_1.fastq
Approx 85% complete for SRR18086277_pass_1.fastq
Approx 90% complete for SRR18086277_pass_1.fastq
Approx 95% complete for SRR18086277_pass_1.fastq

Analysis complete for SRR18086277_pass_1.fastq

Approx 100% complete for SRR18086277_pass_1.fastq
Started analysis of SRR18086277_pass_2.fastq
Approx 5% complete for SRR18086277_pass_2.fastq
Approx 10% complete for SRR18086277_pass_2.fastq
Approx 15% complete for SRR18086277_pass_2.fastq
Approx 20% complete for SRR18086277_pass_2.fastq
Approx 25% complete for SRR18086277_pass_2.fastq
Approx 30% complete for SRR18086277_pass_2.fastq
Approx 35% complete for SRR18086277_pass_2.fastq
Approx 40% complete for SRR18086277_pass_2.fastq
Approx 45% complete for SRR18086277_pass_2.fastq
Approx 50% complete for SRR18086277_pass_2.fastq
Approx 55% complete for SRR18086277_pass_2.fastq
Approx 60% complete for SRR18086277_pass_2.fastq
Approx 65% complete for SRR18086277_pass_2.fastq
Approx 70% complete for SRR18086277_pass_2.fastq
Approx 75% complete for SRR18086277_pass_2.fastq
Approx 80% complete for SRR18086277_pass_2.fastq
Approx 85% complete for SRR18086277_pass_2.fastq
Approx 90% complete for SRR18086277_pass_2.fastq
Approx 95% complete for SRR18086277_pass_2.fastq

Analysis complete for SRR18086277_pass_2.fastq

Approx 100% complete for SRR18086277_pass_2.fastq
Started analysis of SRR18086342_pass_1.fastq
Approx 5% complete for SRR18086342_pass_1.fastq
Approx 10% complete for SRR18086342_pass_1.fastq
Approx 15% complete for SRR18086342_pass_1.fastq
Approx 20% complete for SRR18086342_pass_1.fastq
Approx 25% complete for SRR18086342_pass_1.fastq
Approx 30% complete for SRR18086342_pass_1.fastq
Approx 35% complete for SRR18086342_pass_1.fastq
Approx 40% complete for SRR18086342_pass_1.fastq
Approx 45% complete for SRR18086342_pass_1.fastq
Approx 50% complete for SRR18086342_pass_1.fastq
Approx 55% complete for SRR18086342_pass_1.fastq
Approx 60% complete for SRR18086342_pass_1.fastq
Approx 65% complete for SRR18086342_pass_1.fastq
Approx 70% complete for SRR18086342_pass_1.fastq
Approx 75% complete for SRR18086342_pass_1.fastq
Approx 80% complete for SRR18086342_pass_1.fastq
Approx 85% complete for SRR18086342_pass_1.fastq
Approx 90% complete for SRR18086342_pass_1.fastq
Approx 95% complete for SRR18086342_pass_1.fastq

Analysis complete for SRR18086342_pass_1.fastq

Started analysis of SRR18086342_pass_2.fastq
Approx 5% complete for SRR18086342_pass_2.fastq
Approx 10% complete for SRR18086342_pass_2.fastq
Approx 15% complete for SRR18086342_pass_2.fastq
Approx 20% complete for SRR18086342_pass_2.fastq
Approx 25% complete for SRR18086342_pass_2.fastq
Approx 30% complete for SRR18086342_pass_2.fastq
Approx 35% complete for SRR18086342_pass_2.fastq
Approx 40% complete for SRR18086342_pass_2.fastq
Approx 45% complete for SRR18086342_pass_2.fastq
Approx 50% complete for SRR18086342_pass_2.fastq
Approx 55% complete for SRR18086342_pass_2.fastq
Approx 60% complete for SRR18086342_pass_2.fastq
Approx 65% complete for SRR18086342_pass_2.fastq
Approx 70% complete for SRR18086342_pass_2.fastq
Approx 75% complete for SRR18086342_pass_2.fastq
Approx 80% complete for SRR18086342_pass_2.fastq
Approx 85% complete for SRR18086342_pass_2.fastq
Approx 90% complete for SRR18086342_pass_2.fastq
Approx 95% complete for SRR18086342_pass_2.fastq

Analysis complete for SRR18086342_pass_2.fastq

Started analysis of SRR18086364_pass_1.fastq
Approx 5% complete for SRR18086364_pass_1.fastq
Approx 10% complete for SRR18086364_pass_1.fastq
Approx 15% complete for SRR18086364_pass_1.fastq
Approx 20% complete for SRR18086364_pass_1.fastq
Approx 25% complete for SRR18086364_pass_1.fastq
Approx 30% complete for SRR18086364_pass_1.fastq
Approx 35% complete for SRR18086364_pass_1.fastq
Approx 40% complete for SRR18086364_pass_1.fastq
Approx 45% complete for SRR18086364_pass_1.fastq
Approx 50% complete for SRR18086364_pass_1.fastq
Approx 55% complete for SRR18086364_pass_1.fastq
Approx 60% complete for SRR18086364_pass_1.fastq
Approx 65% complete for SRR18086364_pass_1.fastq
Approx 70% complete for SRR18086364_pass_1.fastq
Approx 75% complete for SRR18086364_pass_1.fastq
Approx 80% complete for SRR18086364_pass_1.fastq
Approx 85% complete for SRR18086364_pass_1.fastq
Approx 90% complete for SRR18086364_pass_1.fastq

Analysis complete for SRR18086364_pass_1.fastq

Approx 95% complete for SRR18086364_pass_1.fastq
Started analysis of SRR18086364_pass_2.fastq
Approx 5% complete for SRR18086364_pass_2.fastq
Approx 10% complete for SRR18086364_pass_2.fastq
Approx 15% complete for SRR18086364_pass_2.fastq
Approx 20% complete for SRR18086364_pass_2.fastq
Approx 25% complete for SRR18086364_pass_2.fastq
Approx 30% complete for SRR18086364_pass_2.fastq
Approx 35% complete for SRR18086364_pass_2.fastq
Approx 40% complete for SRR18086364_pass_2.fastq
Approx 45% complete for SRR18086364_pass_2.fastq
Approx 50% complete for SRR18086364_pass_2.fastq
Approx 55% complete for SRR18086364_pass_2.fastq
Approx 60% complete for SRR18086364_pass_2.fastq
Approx 65% complete for SRR18086364_pass_2.fastq
Approx 70% complete for SRR18086364_pass_2.fastq
Approx 75% complete for SRR18086364_pass_2.fastq
Approx 80% complete for SRR18086364_pass_2.fastq
Approx 85% complete for SRR18086364_pass_2.fastq
Approx 90% complete for SRR18086364_pass_2.fastq

Analysis complete for SRR18086364_pass_2.fastq

Approx 95% complete for SRR18086364_pass_2.fastq

Analysis complete for SRR18086364_pass.fastq

Started analysis of SRR18086364_pass.fastq
Started analysis of SRR18086575_pass_1.fastq
Approx 5% complete for SRR18086575_pass_1.fastq
Approx 10% complete for SRR18086575_pass_1.fastq
Approx 15% complete for SRR18086575_pass_1.fastq
Approx 20% complete for SRR18086575_pass_1.fastq
Approx 25% complete for SRR18086575_pass_1.fastq
Approx 30% complete for SRR18086575_pass_1.fastq
Approx 35% complete for SRR18086575_pass_1.fastq
Approx 40% complete for SRR18086575_pass_1.fastq
Approx 45% complete for SRR18086575_pass_1.fastq
Approx 50% complete for SRR18086575_pass_1.fastq
Approx 55% complete for SRR18086575_pass_1.fastq
Approx 60% complete for SRR18086575_pass_1.fastq
Approx 65% complete for SRR18086575_pass_1.fastq
Approx 70% complete for SRR18086575_pass_1.fastq
Approx 75% complete for SRR18086575_pass_1.fastq
Approx 80% complete for SRR18086575_pass_1.fastq
Approx 85% complete for SRR18086575_pass_1.fastq
Approx 90% complete for SRR18086575_pass_1.fastq
Approx 95% complete for SRR18086575_pass_1.fastq

Analysis complete for SRR18086575_pass_1.fastq

Approx 100% complete for SRR18086575_pass_1.fastq
Started analysis of SRR18086575_pass_2.fastq
Approx 5% complete for SRR18086575_pass_2.fastq
Approx 10% complete for SRR18086575_pass_2.fastq
Approx 15% complete for SRR18086575_pass_2.fastq
Approx 20% complete for SRR18086575_pass_2.fastq
Approx 25% complete for SRR18086575_pass_2.fastq
Approx 30% complete for SRR18086575_pass_2.fastq
Approx 35% complete for SRR18086575_pass_2.fastq
Approx 40% complete for SRR18086575_pass_2.fastq
Approx 45% complete for SRR18086575_pass_2.fastq
Approx 50% complete for SRR18086575_pass_2.fastq
Approx 55% complete for SRR18086575_pass_2.fastq
Approx 60% complete for SRR18086575_pass_2.fastq
Approx 65% complete for SRR18086575_pass_2.fastq
Approx 70% complete for SRR18086575_pass_2.fastq
Approx 75% complete for SRR18086575_pass_2.fastq
Approx 80% complete for SRR18086575_pass_2.fastq
Approx 85% complete for SRR18086575_pass_2.fastq
Approx 90% complete for SRR18086575_pass_2.fastq
Approx 95% complete for SRR18086575_pass_2.fastq

Analysis complete for SRR18086575_pass_2.fastq

Approx 100% complete for SRR18086575_pass_2.fastq
Started analysis of SRR18091316_pass_1.fastq
Approx 5% complete for SRR18091316_pass_1.fastq
Approx 10% complete for SRR18091316_pass_1.fastq
Approx 15% complete for SRR18091316_pass_1.fastq
Approx 20% complete for SRR18091316_pass_1.fastq
Approx 25% complete for SRR18091316_pass_1.fastq
Approx 30% complete for SRR18091316_pass_1.fastq
Approx 35% complete for SRR18091316_pass_1.fastq
Approx 40% complete for SRR18091316_pass_1.fastq
Approx 45% complete for SRR18091316_pass_1.fastq
Approx 50% complete for SRR18091316_pass_1.fastq
Approx 55% complete for SRR18091316_pass_1.fastq
Approx 60% complete for SRR18091316_pass_1.fastq
Approx 65% complete for SRR18091316_pass_1.fastq
Approx 70% complete for SRR18091316_pass_1.fastq
Approx 75% complete for SRR18091316_pass_1.fastq
Approx 80% complete for SRR18091316_pass_1.fastq
Approx 85% complete for SRR18091316_pass_1.fastq
Approx 90% complete for SRR18091316_pass_1.fastq
Approx 95% complete for SRR18091316_pass_1.fastq

Analysis complete for SRR18091316_pass_1.fastq

Approx 100% complete for SRR18091316_pass_1.fastq
Started analysis of SRR18091316_pass_2.fastq
Approx 5% complete for SRR18091316_pass_2.fastq
Approx 10% complete for SRR18091316_pass_2.fastq
Approx 15% complete for SRR18091316_pass_2.fastq
Approx 20% complete for SRR18091316_pass_2.fastq
Approx 25% complete for SRR18091316_pass_2.fastq
Approx 30% complete for SRR18091316_pass_2.fastq
Approx 35% complete for SRR18091316_pass_2.fastq
Approx 40% complete for SRR18091316_pass_2.fastq
Approx 45% complete for SRR18091316_pass_2.fastq
Approx 50% complete for SRR18091316_pass_2.fastq
Approx 55% complete for SRR18091316_pass_2.fastq
Approx 60% complete for SRR18091316_pass_2.fastq
Approx 65% complete for SRR18091316_pass_2.fastq
Approx 70% complete for SRR18091316_pass_2.fastq
Approx 75% complete for SRR18091316_pass_2.fastq
Approx 80% complete for SRR18091316_pass_2.fastq
Approx 85% complete for SRR18091316_pass_2.fastq
Approx 90% complete for SRR18091316_pass_2.fastq
Approx 95% complete for SRR18091316_pass_2.fastq

Analysis complete for SRR18091316_pass_2.fastq

Approx 100% complete for SRR18091316_pass_2.fastq
Started analysis of SRR18091737_pass_1.fastq
Approx 5% complete for SRR18091737_pass_1.fastq
Approx 10% complete for SRR18091737_pass_1.fastq
Approx 15% complete for SRR18091737_pass_1.fastq
Approx 20% complete for SRR18091737_pass_1.fastq
Approx 25% complete for SRR18091737_pass_1.fastq
Approx 30% complete for SRR18091737_pass_1.fastq
Approx 35% complete for SRR18091737_pass_1.fastq
Approx 40% complete for SRR18091737_pass_1.fastq
Approx 45% complete for SRR18091737_pass_1.fastq
Approx 50% complete for SRR18091737_pass_1.fastq
Approx 55% complete for SRR18091737_pass_1.fastq
Approx 60% complete for SRR18091737_pass_1.fastq
Approx 65% complete for SRR18091737_pass_1.fastq
Approx 70% complete for SRR18091737_pass_1.fastq
Approx 75% complete for SRR18091737_pass_1.fastq
Approx 80% complete for SRR18091737_pass_1.fastq
Approx 85% complete for SRR18091737_pass_1.fastq
Approx 90% complete for SRR18091737_pass_1.fastq
Approx 95% complete for SRR18091737_pass_1.fastq

Analysis complete for SRR18091737_pass_1.fastq

Approx 100% complete for SRR18091737_pass_1.fastq
Started analysis of SRR18091737_pass_2.fastq
Approx 5% complete for SRR18091737_pass_2.fastq
Approx 10% complete for SRR18091737_pass_2.fastq
Approx 15% complete for SRR18091737_pass_2.fastq
Approx 20% complete for SRR18091737_pass_2.fastq
Approx 25% complete for SRR18091737_pass_2.fastq
Approx 30% complete for SRR18091737_pass_2.fastq
Approx 35% complete for SRR18091737_pass_2.fastq
Approx 40% complete for SRR18091737_pass_2.fastq
Approx 45% complete for SRR18091737_pass_2.fastq
Approx 50% complete for SRR18091737_pass_2.fastq
Approx 55% complete for SRR18091737_pass_2.fastq
Approx 60% complete for SRR18091737_pass_2.fastq
Approx 65% complete for SRR18091737_pass_2.fastq
Approx 70% complete for SRR18091737_pass_2.fastq
Approx 75% complete for SRR18091737_pass_2.fastq
Approx 80% complete for SRR18091737_pass_2.fastq
Approx 85% complete for SRR18091737_pass_2.fastq
Approx 90% complete for SRR18091737_pass_2.fastq
Approx 95% complete for SRR18091737_pass_2.fastq

Analysis complete for SRR18091737_pass_2.fastq

Approx 100% complete for SRR18091737_pass_2.fastq

#source ~/.bashrc
module load multiqc/1.7.0
mkdir -p multiqc
multiqc ./fastqc -o ./multiqc

[WARNING]         multiqc : MultiQC Version v1.13 now available!
[INFO   ]         multiqc : This is MultiQC v1.7
[INFO   ]         multiqc : Template    : default
[INFO   ]         multiqc : Searching './fastqc'
[INFO   ]          fastqc : Found 17 reports
[INFO   ]         multiqc : Compressing plot data
[INFO   ]         multiqc : Report      : multiqc/multiqc_report.html
[INFO   ]         multiqc : Data        : multiqc/multiqc_data
[INFO   ]         multiqc : MultiQC complete

Check the quality of your results by navigating to the file folder for multiqc and clicking on the html files that were generated.

Okay, so the data is not the most perfect, but we just need a set that can make it through the analysis.

Preprocessing Reads

Metagenomics has benefited from k-mer based approaches.

BCBio tools are useful for cleaining metagenomics files, because many simple functions can be called and sped up using kmers.

For example, let’s remove known adapter and other contamination using the NCBI database called “UniVec”

cd ~/Metagenomics2022
#pwd

#grep "TruSeq" ./contaminants/UniVec.fasta

grep "PCR Primer Index 9" ./contaminants/UniVec.fasta

>gnl|uv|NGB00802.1:1-63 Illumina TruSeq Small RNA Sample Prep Kit RNA PCR Primer Index 9 (RPI9) (Oligonucleotide sequence copyright 2007-2012 Illumina, Inc. All rights reserved.)
>gnl|uv|NGB00372.1:1-43 Illumina Multiplexing PCR Primer Index 9 (Oligonucleotide sequence copyright 2007-2012 Illumina, Inc. All rights reserved.)

Remove contamination from sequences with BCBio tools

Are these repeats biological (part of the genome) or technical (primer amplification, underabundance of sample or lack of PhiX spike-in)? What does the raw sequencing data show us?

Here is an example of how bbduk.sh can be used to simultaneously remove adapters, phiX spike-in and overrepresented sequences from our custom file.

#source ~/.bashrc
module load java/1.8.0_60
module load /cluster/tufts/bio/tools/module_files/bcbio/1.1.5
mkdir -p decon


bbduk.sh -Xmx10g -tbo  -tpe \
in1=data/SRR18085691_pass_1.fastq \
in2=data/SRR18085691_pass_2.fastq \
out1=decon/SRR18085691_decon_1.fastq \
out2=decon/SRR18085691_decon_2.fastq \
outs=decon/SRR18085691_decon_sing.fastq \
ref=contaminants/UniVec.fasta \
ref=contaminants/overrep.fasta \
ref=phix \
ktrim=r k=21 rcomp=t mink=11 hdist=2 minlen=200 \

java -Djava.library.path=/cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/jni/ -ea -Xmx10g -Xms10g -cp /cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/current/ jgi.BBDuk -Xmx10g -tbo -tpe in1=data/SRR18085691_pass_1.fastq in2=data/SRR18085691_pass_2.fastq out1=decon/SRR18085691_decon_1.fastq out2=decon/SRR18085691_decon_2.fastq outs=decon/SRR18085691_decon_sing.fastq ref=contaminants/UniVec.fasta ref=contaminants/overrep.fasta ref=phix ktrim=r k=21 rcomp=t mink=11 hdist=2 minlen=200
Executing jgi.BBDuk [-Xmx10g, tbo, tpe, in1=data/SRR18085691_pass_1.fastq, in2=data/SRR18085691_pass_2.fastq, out1=decon/SRR18085691_decon_1.fastq, out2=decon/SRR18085691_decon_2.fastq, outs=decon/SRR18085691_decon_sing.fastq, ref=contaminants/UniVec.fasta, ref=contaminants/overrep.fasta, ref=phix, ktrim=r, k=21, rcomp=t, mink=11, hdist=2, minlen=200]
Version 38.46

maskMiddle was disabled because useShortKmers=true
0.036 seconds.
Initial:
Memory: max=10290m, total=10290m, free=9968m, used=322m

Added 10506850 kmers; time:     5.863 seconds.
Memory: max=10290m, total=10290m, free=9233m, used=1057m

Input is being processed as paired
Started output streams: 0.066 seconds.
Processing time:        81.613 seconds.

Input:                      2000000 reads       502000000 bases.
KTrimmed:                   10970 reads (0.55%)     1959602 bases (0.39%)
Trimmed by overlap:         112718 reads (5.64%)    6611950 bases (1.32%)
Total Removed:              57438 reads (2.87%)     15942726 bases (3.18%)
Result:                     1942562 reads (97.13%)  486057274 bases (96.82%)

Time:                           87.544 seconds.
Reads Processed:       2000k    22.85k reads/sec
Bases Processed:        502m    5.73m bases/sec

While the next chunk is running, let’s talk about kmer-based approaches to filtering.

BBDUK = Decontamination Using Kmers

• -tbo trims overlapping bases in paired end reads
• -tpe trims both paired end reads to the same length
• -in1,-in2,out1,out2,outs means input forward and reverse reads and output forward and reverse plus orphaned reads
• ref can be any file in fasta format or baked-in options such as phiX and adapters
• ref can be replaced with literal=ACTGGT,TTTGGTG with any list of strings
• ktrim=r means to trim from 3-prime end
• k refers to the kmer size that is chosen
• mink allows a shorter kmer size at the end of the sequence
• hdist refers to hamming distance - number of mismatches allowed in a kmer
• minlen indicates the minimum length of the sequence to keep after trimming

#source ~/.bashrc
module load java/1.8.0_60
module load /cluster/tufts/bio/tools/module_files/bcbio/1.1.5
mkdir -p decon


ls data/*_pass_1.fastq | while read line ; do \

bbduk.sh -Xmx10g -tbo -tpe \
in1=$line \
in2=data/$(basename $line _pass_1.fastq)_pass_2.fastq \
out1=decon/$(basename $line _pass_1.fastq)_decon_1.fastq \
out2=decon/$(basename $line _pass_1.fastq)_decon_2.fastq\ \
outs=decon/$(basename $line _pass_1.fastq)_decon_sing.fastq \
ref=contaminants/UniVec.fasta \
ref=contaminants/overrep.fasta \
ref=phix \
ktrim=r k=21 rcomp=t mink=11 hdist=2 minlen=200; \
done

java -Djava.library.path=/cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/jni/ -ea -Xmx10g -Xms10g -cp /cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/current/ jgi.BBDuk -Xmx10g -tbo -tpe in1=data/SRR18085691_pass_1.fastq in2=data/SRR18085691_pass_2.fastq out1=decon/SRR18085691_decon_1.fastq out2=decon/SRR18085691_decon_2.fastq outs=decon/SRR18085691_decon_sing.fastq ref=contaminants/UniVec.fasta ref=contaminants/overrep.fasta ref=phix ktrim=r k=21 rcomp=t mink=11 hdist=2 minlen=200
Executing jgi.BBDuk [-Xmx10g, tbo, tpe, in1=data/SRR18085691_pass_1.fastq, in2=data/SRR18085691_pass_2.fastq, out1=decon/SRR18085691_decon_1.fastq, out2=decon/SRR18085691_decon_2.fastq, outs=decon/SRR18085691_decon_sing.fastq, ref=contaminants/UniVec.fasta, ref=contaminants/overrep.fasta, ref=phix, ktrim=r, k=21, rcomp=t, mink=11, hdist=2, minlen=200]
Version 38.46

maskMiddle was disabled because useShortKmers=true
0.037 seconds.
Initial:
Memory: max=10290m, total=10290m, free=9859m, used=431m

Added 10506850 kmers; time:     5.765 seconds.
Memory: max=10290m, total=10290m, free=9124m, used=1166m

Input is being processed as paired
Started output streams: 0.374 seconds.
Processing time:        83.866 seconds.

Input:                      2000000 reads       502000000 bases.
KTrimmed:                   10970 reads (0.55%)     1959602 bases (0.39%)
Trimmed by overlap:         112718 reads (5.64%)    6611950 bases (1.32%)
Total Removed:              57438 reads (2.87%)     15942726 bases (3.18%)
Result:                     1942562 reads (97.13%)  486057274 bases (96.82%)

Time:                           90.007 seconds.
Reads Processed:       2000k    22.22k reads/sec
Bases Processed:        502m    5.58m bases/sec
java -Djava.library.path=/cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/jni/ -ea -Xmx10g -Xms10g -cp /cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/current/ jgi.BBDuk -Xmx10g -tbo -tpe in1=data/SRR18085708_pass_1.fastq in2=data/SRR18085708_pass_2.fastq out1=decon/SRR18085708_decon_1.fastq out2=decon/SRR18085708_decon_2.fastq outs=decon/SRR18085708_decon_sing.fastq ref=contaminants/UniVec.fasta ref=contaminants/overrep.fasta ref=phix ktrim=r k=21 rcomp=t mink=11 hdist=2 minlen=200
Executing jgi.BBDuk [-Xmx10g, tbo, tpe, in1=data/SRR18085708_pass_1.fastq, in2=data/SRR18085708_pass_2.fastq, out1=decon/SRR18085708_decon_1.fastq, out2=decon/SRR18085708_decon_2.fastq, outs=decon/SRR18085708_decon_sing.fastq, ref=contaminants/UniVec.fasta, ref=contaminants/overrep.fasta, ref=phix, ktrim=r, k=21, rcomp=t, mink=11, hdist=2, minlen=200]
Version 38.46

maskMiddle was disabled because useShortKmers=true
0.030 seconds.
Initial:
Memory: max=10290m, total=10290m, free=9914m, used=376m

Added 10506850 kmers; time:     6.912 seconds.
Memory: max=10290m, total=10290m, free=9179m, used=1111m

Input is being processed as paired
Started output streams: 0.058 seconds.
Processing time:        70.165 seconds.

Input:                      2000000 reads       502000000 bases.
KTrimmed:                   2840 reads (0.14%)  360174 bases (0.07%)
Trimmed by overlap:         195124 reads (9.76%)    25710564 bases (5.12%)
Total Removed:              161564 reads (8.08%)    41435452 bases (8.25%)
Result:                     1838436 reads (91.92%)  460564548 bases (91.75%)

Time:                           77.137 seconds.
Reads Processed:       2000k    25.93k reads/sec
Bases Processed:        502m    6.51m bases/sec
java -Djava.library.path=/cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/jni/ -ea -Xmx10g -Xms10g -cp /cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/current/ jgi.BBDuk -Xmx10g -tbo -tpe in1=data/SRR18086277_pass_1.fastq in2=data/SRR18086277_pass_2.fastq out1=decon/SRR18086277_decon_1.fastq out2=decon/SRR18086277_decon_2.fastq outs=decon/SRR18086277_decon_sing.fastq ref=contaminants/UniVec.fasta ref=contaminants/overrep.fasta ref=phix ktrim=r k=21 rcomp=t mink=11 hdist=2 minlen=200
Executing jgi.BBDuk [-Xmx10g, tbo, tpe, in1=data/SRR18086277_pass_1.fastq, in2=data/SRR18086277_pass_2.fastq, out1=decon/SRR18086277_decon_1.fastq, out2=decon/SRR18086277_decon_2.fastq, outs=decon/SRR18086277_decon_sing.fastq, ref=contaminants/UniVec.fasta, ref=contaminants/overrep.fasta, ref=phix, ktrim=r, k=21, rcomp=t, mink=11, hdist=2, minlen=200]
Version 38.46

maskMiddle was disabled because useShortKmers=true
0.039 seconds.
Initial:
Memory: max=10290m, total=10290m, free=9860m, used=430m

Added 10506850 kmers; time:     6.574 seconds.
Memory: max=10290m, total=10290m, free=9126m, used=1164m

Input is being processed as paired
Started output streams: 0.040 seconds.
Processing time:        79.791 seconds.

Input:                      2000000 reads       502000000 bases.
KTrimmed:                   918 reads (0.05%)   198378 bases (0.04%)
Trimmed by overlap:         141132 reads (7.06%)    22658800 bases (4.51%)
Total Removed:              134148 reads (6.71%)    33887548 bases (6.75%)
Result:                     1865852 reads (93.29%)  468112452 bases (93.25%)

Time:                           86.407 seconds.
Reads Processed:       2000k    23.15k reads/sec
Bases Processed:        502m    5.81m bases/sec
java -Djava.library.path=/cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/jni/ -ea -Xmx10g -Xms10g -cp /cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/current/ jgi.BBDuk -Xmx10g -tbo -tpe in1=data/SRR18086342_pass_1.fastq in2=data/SRR18086342_pass_2.fastq out1=decon/SRR18086342_decon_1.fastq out2=decon/SRR18086342_decon_2.fastq outs=decon/SRR18086342_decon_sing.fastq ref=contaminants/UniVec.fasta ref=contaminants/overrep.fasta ref=phix ktrim=r k=21 rcomp=t mink=11 hdist=2 minlen=200
Executing jgi.BBDuk [-Xmx10g, tbo, tpe, in1=data/SRR18086342_pass_1.fastq, in2=data/SRR18086342_pass_2.fastq, out1=decon/SRR18086342_decon_1.fastq, out2=decon/SRR18086342_decon_2.fastq, outs=decon/SRR18086342_decon_sing.fastq, ref=contaminants/UniVec.fasta, ref=contaminants/overrep.fasta, ref=phix, ktrim=r, k=21, rcomp=t, mink=11, hdist=2, minlen=200]
Version 38.46

maskMiddle was disabled because useShortKmers=true
0.027 seconds.
Initial:
Memory: max=10290m, total=10290m, free=9966m, used=324m

Added 10506850 kmers; time:     7.036 seconds.
Memory: max=10290m, total=10290m, free=9232m, used=1058m

Input is being processed as paired
Started output streams: 0.036 seconds.
Processing time:        55.800 seconds.

Input:                      1334530 reads       334967030 bases.
KTrimmed:                   7194 reads (0.54%)  1193470 bases (0.36%)
Trimmed by overlap:         95964 reads (7.19%)     7114482 bases (2.12%)
Total Removed:              58776 reads (4.40%)     15820388 bases (4.72%)
Result:                     1275754 reads (95.60%)  319146642 bases (95.28%)

Time:                           62.875 seconds.
Reads Processed:       1334k    21.23k reads/sec
Bases Processed:        334m    5.33m bases/sec
java -Djava.library.path=/cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/jni/ -ea -Xmx10g -Xms10g -cp /cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/current/ jgi.BBDuk -Xmx10g -tbo -tpe in1=data/SRR18086364_pass_1.fastq in2=data/SRR18086364_pass_2.fastq out1=decon/SRR18086364_decon_1.fastq out2=decon/SRR18086364_decon_2.fastq outs=decon/SRR18086364_decon_sing.fastq ref=contaminants/UniVec.fasta ref=contaminants/overrep.fasta ref=phix ktrim=r k=21 rcomp=t mink=11 hdist=2 minlen=200
Executing jgi.BBDuk [-Xmx10g, tbo, tpe, in1=data/SRR18086364_pass_1.fastq, in2=data/SRR18086364_pass_2.fastq, out1=decon/SRR18086364_decon_1.fastq, out2=decon/SRR18086364_decon_2.fastq, outs=decon/SRR18086364_decon_sing.fastq, ref=contaminants/UniVec.fasta, ref=contaminants/overrep.fasta, ref=phix, ktrim=r, k=21, rcomp=t, mink=11, hdist=2, minlen=200]
Version 38.46

maskMiddle was disabled because useShortKmers=true
0.031 seconds.
Initial:
Memory: max=10290m, total=10290m, free=9913m, used=377m

Added 10506850 kmers; time:     6.257 seconds.
Memory: max=10290m, total=10290m, free=9178m, used=1112m

Input is being processed as paired
Started output streams: 0.156 seconds.
Processing time:        5.335 seconds.

Input:                      59630 reads         14875247 bases.
KTrimmed:                   138 reads (0.23%)   20946 bases (0.14%)
Trimmed by overlap:         10 reads (0.02%)    283 bases (0.00%)
Total Removed:              84 reads (0.14%)    21591 bases (0.15%)
Result:                     59546 reads (99.86%)    14853656 bases (99.85%)

Time:                           11.750 seconds.
Reads Processed:       59630    5.08k reads/sec
Bases Processed:      14875k    1.27m bases/sec
java -Djava.library.path=/cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/jni/ -ea -Xmx10g -Xms10g -cp /cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/current/ jgi.BBDuk -Xmx10g -tbo -tpe in1=data/SRR18086575_pass_1.fastq in2=data/SRR18086575_pass_2.fastq out1=decon/SRR18086575_decon_1.fastq out2=decon/SRR18086575_decon_2.fastq outs=decon/SRR18086575_decon_sing.fastq ref=contaminants/UniVec.fasta ref=contaminants/overrep.fasta ref=phix ktrim=r k=21 rcomp=t mink=11 hdist=2 minlen=200
Executing jgi.BBDuk [-Xmx10g, tbo, tpe, in1=data/SRR18086575_pass_1.fastq, in2=data/SRR18086575_pass_2.fastq, out1=decon/SRR18086575_decon_1.fastq, out2=decon/SRR18086575_decon_2.fastq, outs=decon/SRR18086575_decon_sing.fastq, ref=contaminants/UniVec.fasta, ref=contaminants/overrep.fasta, ref=phix, ktrim=r, k=21, rcomp=t, mink=11, hdist=2, minlen=200]
Version 38.46

maskMiddle was disabled because useShortKmers=true
0.021 seconds.
Initial:
Memory: max=10290m, total=10290m, free=9859m, used=431m

Added 10506850 kmers; time:     6.260 seconds.
Memory: max=10290m, total=10290m, free=9124m, used=1166m

Input is being processed as paired
Started output streams: 0.035 seconds.
Processing time:        92.165 seconds.

Input:                      2000000 reads       502000000 bases.
KTrimmed:                   3002 reads (0.15%)  460100 bases (0.09%)
Trimmed by overlap:         439534 reads (21.98%)   62037664 bases (12.36%)
Total Removed:              381196 reads (19.06%)   97227296 bases (19.37%)
Result:                     1618804 reads (80.94%)  404772704 bases (80.63%)

Time:                           98.463 seconds.
Reads Processed:       2000k    20.31k reads/sec
Bases Processed:        502m    5.10m bases/sec
java -Djava.library.path=/cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/jni/ -ea -Xmx10g -Xms10g -cp /cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/current/ jgi.BBDuk -Xmx10g -tbo -tpe in1=data/SRR18091316_pass_1.fastq in2=data/SRR18091316_pass_2.fastq out1=decon/SRR18091316_decon_1.fastq out2=decon/SRR18091316_decon_2.fastq outs=decon/SRR18091316_decon_sing.fastq ref=contaminants/UniVec.fasta ref=contaminants/overrep.fasta ref=phix ktrim=r k=21 rcomp=t mink=11 hdist=2 minlen=200
Executing jgi.BBDuk [-Xmx10g, tbo, tpe, in1=data/SRR18091316_pass_1.fastq, in2=data/SRR18091316_pass_2.fastq, out1=decon/SRR18091316_decon_1.fastq, out2=decon/SRR18091316_decon_2.fastq, outs=decon/SRR18091316_decon_sing.fastq, ref=contaminants/UniVec.fasta, ref=contaminants/overrep.fasta, ref=phix, ktrim=r, k=21, rcomp=t, mink=11, hdist=2, minlen=200]
Version 38.46

maskMiddle was disabled because useShortKmers=true
0.031 seconds.
Initial:
Memory: max=10290m, total=10290m, free=9860m, used=430m

Added 10506850 kmers; time:     6.352 seconds.
Memory: max=10290m, total=10290m, free=9126m, used=1164m

Input is being processed as paired
Started output streams: 0.098 seconds.
Processing time:        88.557 seconds.

Input:                      2000000 reads       502000000 bases.
KTrimmed:                   11350 reads (0.57%)     2233876 bases (0.44%)
Trimmed by overlap:         141514 reads (7.08%)    9911700 bases (1.97%)
Total Removed:              85398 reads (4.27%)     22989670 bases (4.58%)
Result:                     1914602 reads (95.73%)  479010330 bases (95.42%)

Time:                           95.010 seconds.
Reads Processed:       2000k    21.05k reads/sec
Bases Processed:        502m    5.28m bases/sec
java -Djava.library.path=/cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/jni/ -ea -Xmx10g -Xms10g -cp /cluster/tufts/bio/tools/bcbio/1.1.5/anaconda/opt/bbmap-38.46-0/current/ jgi.BBDuk -Xmx10g -tbo -tpe in1=data/SRR18091737_pass_1.fastq in2=data/SRR18091737_pass_2.fastq out1=decon/SRR18091737_decon_1.fastq out2=decon/SRR18091737_decon_2.fastq outs=decon/SRR18091737_decon_sing.fastq ref=contaminants/UniVec.fasta ref=contaminants/overrep.fasta ref=phix ktrim=r k=21 rcomp=t mink=11 hdist=2 minlen=200
Executing jgi.BBDuk [-Xmx10g, tbo, tpe, in1=data/SRR18091737_pass_1.fastq, in2=data/SRR18091737_pass_2.fastq, out1=decon/SRR18091737_decon_1.fastq, out2=decon/SRR18091737_decon_2.fastq, outs=decon/SRR18091737_decon_sing.fastq, ref=contaminants/UniVec.fasta, ref=contaminants/overrep.fasta, ref=phix, ktrim=r, k=21, rcomp=t, mink=11, hdist=2, minlen=200]
Version 38.46

maskMiddle was disabled because useShortKmers=true
0.026 seconds.
Initial:
Memory: max=10290m, total=10290m, free=9859m, used=431m

Added 10506850 kmers; time:     6.520 seconds.
Memory: max=10290m, total=10290m, free=9124m, used=1166m

Input is being processed as paired
Started output streams: 0.074 seconds.
Processing time:        75.760 seconds.

Input:                      2000000 reads       502000000 bases.
KTrimmed:                   768 reads (0.04%)   177766 bases (0.04%)
Trimmed by overlap:         118910 reads (5.95%)    19108262 bases (3.81%)
Total Removed:              114844 reads (5.74%)    28963734 bases (5.77%)
Result:                     1885156 reads (94.26%)  473036266 bases (94.23%)

Time:                           82.357 seconds.
Reads Processed:       2000k    24.28k reads/sec
Bases Processed:        502m    6.10m bases/sec

#source ~/.bashrc
module load fastqc
mkdir -p fastqc_decon
fastqc ./decon/*decon*.fastq -o fastqc_decon

#source ~/.bashrc
module load multiqc/1.7.0
mkdir -p multiqc_decon
multiqc ./fastqc_decon -o ./multiqc_decon

Question: Did we get rid of all the adapter sequences? What would your next step be?

Question: Is there a particular tick pool that we may want to remove before assembly?

To be removed before assembly:

• Ignore single reads
• Ignore low count reads (SRR18086364)
• Ignore low quality reads (SRR18091737,SRR18086277)

mkdir -p gooddata

echo SRR18085691 > gooddata.txt
echo SRR18085708 >> gooddata.txt
echo SRR18086342 >> gooddata.txt
echo SRR18086575 >> gooddata.txt
echo SRR18091316 >> gooddata.txt


cat gooddata.txt | while read $line; do cp decon/$line*1.fastq gooddata; cp decon/$line*2.fastq gooddata; done

This concludes Part I, please come back after the break for Part II - using a wrapper script to bin your data.

---

Part II

Running an Assembly on the Cleaned Data

This command takes approximately 20 minutes, so the commands used to make the assembly are commented out.

#/cluster/tufts/bio/tools/megahit/1.2.9/bin/megahit  -m 16000000 -t 32 --read ./gooddata/* \
#--k-list 21,41,61,81,99 \
#--no-mercy \
#--min-count 2 \
#-o megahit_out/

Megahit: de Bruijn graph Assembly

Megahit workflow

What is a deBruijn graph?

It is an assembly algorithm that splits your sequence into k-mers and tries to build the longest contigs along the shortest parsimonious path.

Megahit workflow

To save time in the workshop, the assembly has been done using the commented out command above and placed into your workshop folder called “reports”.

Checking the quality of our Assembly

Quast is a commonly used tool to assess assemblies. It will provide a report that is either in a text file or an html.

Let’s run the following command and then look for the quast report to see if we have a “good” assembly.

/cluster/tufts/bio/tools/quast/5.2.0/quast.py -o quast_assembly ./reports/final.contigs.fa

/cluster/tufts/bio/tools/quast/5.2.0/quast.py -o quast_assembly ./reports/final.contigs.fa

Version: 5.2.0, 3d87c606

System information:
  OS: Linux-3.10.0-1127.el7.x86_64-x86_64-with-redhat-7.8-Maipo (linux_64)
  Python version: 2.7.5
  CPUs number: 72

Started: 2022-11-09 11:19:20

Logging to /cluster/home/arhode05/Metagenomics2022/quast_assembly/quast.log
NOTICE: Maximum number of threads is set to 18 (use --threads option to set it manually)

CWD: /cluster/home/arhode05/Metagenomics2022
Main parameters: 
  MODE: default, threads: 18, min contig length: 500, min alignment length: 65, min alignment IDY: 95.0, \
  ambiguity: one, min local misassembly length: 200, min extensive misassembly length: 1000

WARNING: Can't draw plots: python-matplotlib is missing or corrupted.

Contigs:
  Pre-processing...
  ./reports/final.contigs.fa ==> final.contigs

2022-11-09 11:19:22
Running Basic statistics processor...
  Contig files: 
    final.contigs
  Calculating N50 and L50...
    final.contigs, N50 = 876, L50 = 7401, auN = 2787.5, Total length = 21673500, GC % = 53.33, # N's per 100 kbp =  0.00
Done.

NOTICE: Genes are not predicted by default. Use --gene-finding or --glimmer option to enable it.

2022-11-09 11:19:23
Creating large visual summaries...
This may take a while: press Ctrl-C to skip this step..
  1 of 1: Creating Icarus viewers...
Done

2022-11-09 11:19:24
RESULTS:
  Text versions of total report are saved to /cluster/home/arhode05/Metagenomics2022/quast_assembly/report.txt, report.tsv, and report.tex
  Text versions of transposed total report are saved to /cluster/home/arhode05/Metagenomics2022/quast_assembly/transposed_report.txt, transposed_report.tsv, and transposed_report.tex
  HTML version (interactive tables and plots) is saved to /cluster/home/arhode05/Metagenomics2022/quast_assembly/report.html
  Icarus (contig browser) is saved to /cluster/home/arhode05/Metagenomics2022/quast_assembly/icarus.html
  Log is saved to /cluster/home/arhode05/Metagenomics2022/quast_assembly/quast.log

Finished: 2022-11-09 11:19:24
Elapsed time: 0:00:04.007177
NOTICEs: 2; WARNINGs: 1; non-fatal ERRORs: 0

Thank you for using QUAST!

Using Kraken to check for Metagenome Composition

At this point, we may be interested in whether our assembly makes sense for our set of contigs.

Kraken2 is the newest version of Kraken, a taxonomic classification system using exact k-mer matches. This k-mer based approach has very fast classification speeds with high accuracy.

This approach differs from homology based approaches that try to match sequences to each other and score them based on the number of mismatches, deletions and inserts. Kraken uses the entire kmer composition of the contig to match it to a database where reference sequences have been broken down into kmer “hashes” or a database of what the kmers look like for a particular organism.

/cluster/tufts/bio/tools/conda_envs/kraken/2.1.2/bin/kraken2 --use-names --threads 4 --db /cluster/tufts/bio/tools/training/metagenomics/kraken_virus --report kraken.report.txt --quick reports/final.contigs.fa --classified-out classified_sequences.tsv > sequences.kraken

Loading database information... done.
48954 sequences (31.42 Mbp) processed in 0.473s (6215.0 Kseq/m, 3988.65 Mbp/m).
  44345 sequences classified (90.59%)
  4609 sequences unclassified (9.41%)

head classified_sequences.tsv

>k141_13788 flag=1 multi=2.0000 len=335 kraken:taxid|283
TGACGGCGAAGAGGGCGAGGTCAACCGCCTGCAGTTTGCACCGCCCTTTGGCGAGAAGGAGACGATGCGCATCGTGCTGCCCCAGGGCTTCAAGGACGAATCGGACCGGCCGCTGGCCAATGCCAACCTGTTCCCGCTGCAGGTGCAAAGCGGCCTGATGCCGCCGCTGGCCAAGTTTGCCGCCTCGCCCTTTGGCGTGGTCGAGCGCTTTGCCGAAGGCGACAAGGGCCCGGCGATGTTCCCGGTGACCCTGCGCAATGTGGAGCAGGACCTGCAGATCCAGGGCTTGTCGATCGAGGACAAGGGCACCGTCTCGACCCTGCAGCTGCAAAGCG
>k141_36756 flag=1 multi=4.0000 len=1049 kraken:taxid|568
CCTCATGGAGATTTACCAAATGGCTAAGGCAAAACTTAACCCGTTAGCAACTTTATTGCGCGGGGGGCTTTGCGTACCGTTATTCATCGCAGCAGCACCAGTGATGGCCGCGGAAGGGGACGACGTGATGGTTGTCACCGCTTCAGCAACAGAACAAAACCTGAAAGATGCACCTGCCAGTATTAGCGTGATTACGCGTGAAGACTTGGCTAAAAAGCCGATTCAGAACCTGAAGGACATTCTGCAAGAAGTCCCAGGCGTACAGCTGACTAACGAAGGTGATAACCGCAAGGGCGTGAGTATCCGTGGTTTGGACAGTAGCTATACGCTGATTTTGGTGGATGGCAAACGCGTCAGCTCACGTACGGCGGTGTTCCGTCATAACGACTATGACCTAAGCTGGGTGCCCGCCGAGGCGATTGAGCGTATCGAAGTGGTGCGCGGCCCAATGTCTTCGCTTTACGGCTCTGATGCTTTAGGCGGCGTAGTGAACATTATTACCCGTAAAGTGGGTAAAGAGTGGCACGGTTCCCTGAGTGCCGATAGCACTATTCAGGAACACCGCGATCGCGGCGATAGCTATAACGGCAGCTTCTATACCAGTGGCCCGTTGATTGATGATGTGCTCGGAGTAAAAGTCTTTGGCAATTTGGGCAAGCGCGAAAAAGATGATGCCCAGAAATCGACGACCAGCAAATCGGGCGAATCGGGGCGCATAGAAGGGTACACCACGCGTGATGCTAACGTTGAGTTCGCATGGACGCCAGAAGCGAATCAGGACATGACCTTTGGCTACGGTTTTGACCGTCAGGATCGTGATTCTGATTCTTTGGATAAAAACCGTCTGGAACGTCAAAACTATTCAATTGGTCACAACGGGCGTTGGGATTTAGCGAACACCGAAATTCGCGCCTATGGTGAAAAGATTGATAACTACAACGAAAATCGTATCACGGCGAAAAATAACGCCATCGACGGGAAAGTTGTTGTGCCCCTCGGTGAGCTGAATCAGCTATTCACACTGGGTGGTGAATACCGTAACGATAAAC
>k141_36757 flag=1 multi=2.0000 len=346 kraken:taxid|1907578
TCTCATGGAAGCGTTCCAGTAAGGCTTCAGCGGCCAGCAGCGAGAAACCGCGCGCTTCACGGTTGCCGGTCAGACGCCAGGCCGCCACGCACGACAGGGTATACAGAAAGCCCAAATCGTGAGTATTGGTGTCGTTGCGCCCGGCGATGCGCAGGCCAAACGAGCGCACATTTTTCTCGGCAAGCTGGCGGAATTTTTCGTCGCCGCTCATCTCCCACGCCAGCCACAGCTGCCCGGTCCAGAAGCTGGTGGTCCACTCGACGTTGTCCGTCAGCGGGTAGAACCCATCCTCACAGGTTTCCGCCGGGAATTTTTCACCGAACTCCGTTAAATGACGGCTAATCAG
>k141_13790 flag=1 multi=2.0000 len=331 kraken:taxid|286
CGCTTTGCACTTGCACGTCGATGGCGTTTTCCAAGCCGTCGGGGAACAGGTAGTGGTGGTCGGCCGCCGTGGTGCAGCCCGACAGCAACAATTCGGCCAGCGCGACTTTAGTGGCGAGGGCGAGTTTTTCCGGGGTAAGCCGGGCCCAGACCGGGTACAGGGTTTTCAGCCACGGGAACAACGGTTGGTTGACCACCGGCGCCCAGGCGCGGGTCAGGGTTTGATAGAAGTGGTGATGGGTGTTGATCAGGCCCGGCAGGATTACGTGCTGGCGAGCGTCGAATATGTTTTCGCAAGGTACCGAGGGTTCGCGTCCCCAGCCCAGCACTTC
>k141_13791 flag=1 multi=3.0000 len=771 kraken:taxid|106649
GCAATCTTGTTTATATTTCATGGCATTTTTCTGTTGTTTTTATTGAGACCTCAATCAAAACGCAATTTAGCTCAAAATGCAGATGAATATTTGATTGTATATGCAAGTCAGTCAGGGCAGACAGAAAGTTTTGCTTTGCAAATAGCCCAGCAATTATCTGCTTCAGGTTATTTGGTTCAGTGCTTGAGTATTGAACAAGTGAATGCTGCTTTATTCAAACATGGGACAAAGATTTTATGGATCGTCAGTACTTATGGCGAGGGGGATGCACCAGATTCGGCACGTCATTTTGTTGGGAATGTAATGACAAAGTCATTTGATTTAAATCACATTGAGTATGCTGTGTTGGCTTTTGGTGATTCGCATTATGCAAATTTCTGTAATTTTGGTAGACAATTGGATTTATGGTTAACTGAAAACCATGCACAAGCTTTATTTCCCATGATTTGTGTTGATCAACTAAACCCTACAGGTCTGACGGAATGGCAAACTCAACTTAACCTAGTTTTGCCGCAATTTTCCGTTCAATCTGTTCAATTGAAACAATATCCTTTTCATTCCGTGGTTTTATCTGAACGCCATTTGTTAAATCAAGGCAGTTTAGGGCATCCAATGTATCATCTTGTGTTTACCCAAGTTGAAGCATTATGTTGGAAATCTGGTGATATTCTTGAAATTCAATGTGCCAATACAGATTCACAAATCAAGGCATTTTTAGATCATTTTAATCTTCAAAACTCCAGTGTGTTTCAGGCTAATCATCCAATAGCG

We can open up the sequence report inside RStudio OnDemand and scroll through the findings.

This can be a very long list of species names, so it may be more helpful to use an interactive visualization tool to look at our data.

Krona is a very useful tool that can take in a variety of formats. In this case, we are providing our kraken.report to the visualization. Metabarcoding data assigned to a taxonomy can also be used in this case.

/cluster/tufts/bio/tools/conda_envs/kraken/2.1.2/bin/ktImportTaxonomy -t 5 -m 3 -o krona.html kraken.report.txt

   [ WARNING ]  Score column already in use; not reading scores.
   [ WARNING ]  The following taxonomy IDs were not found in the local
                database and were set to root (if they were recently added to
                NCBI, use updateTaxonomy.sh to update the local database):
                28883 541000 2705459 143813 2546450 741091

Loading taxonomy...
Importing kraken.report.txt...
Writing krona.html...

The output, when opened in a web browser, allows the user to interact with the data and explore different levels of the data.

Krona output

What can we do with Metagenomic Approaches

Case Study: Describing a Complete Cyanobacterial Genome and its Symbionts from an Extremely Arid Environment

The Atacama Desert - Are we on Mars?

Atacama Desert

So we are probably running out of time to continue working hands on.

Let’s take some time to discuss what the next steps might look like.

This discussion will center around a recent project from Dr. Pearson’s lab, conducted by Neveda Naz, Ph.D.

Neveda builds little robots and sends them to outerspace. (Neveda if you are here, feel free to comment!)

Please open up this pdf here to find the discussion.

Helpful Resources

Also, for future reading, I have placed two useful .pdf’s in our public github pages that summarize the steps we are discussing in further detail.

Essentials in Metagenomics I Essentials in Metagenomics II

This markdown and associated files will be placed in our Github and mailed out to the class shortly.

We will also send out a video of this workshop.

Please Provide Feedback

Jason and I hope to continue developing this practical material to provide more hands-on experience in the future, so any feedback on the tools you would like to explore is appreciated.